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We present calculations of the electronic structure of various atoms and molecules in strong 
magnetic fields ranging from B = 10 12 G to 2 x 10 15 G, appropriate for radio pulsars and magnetars. 
For these field strengths, the magnetic forces on the electrons dominate over the Coulomb forces, and 
to a good approximation the electrons are confined to the ground Landau level. Our calculations are 
based on the density functional theory, and use a local magnetic exchange-correlation function which 
is tested to be reliable in the strong field regime. Numerical results of the ground-state energies are 
given for Hjv (up to N = 10), Ue N (up to N = 8), Cjv (up to TV = 5), and Fe N (up to N = 3), 
as well as for various ionized atoms. Fitting formulae for the 73-dependence of the energies are also 
given. In general, as N increases, the binding energy per atom in a molecule, \En\/N, increases and 
approaches a constant value. For all the field strengths considered in this paper, hydrogen, helium, 
and carbon molecules are found to be bound relative to individual atoms (although for B less than a 
few xlO 12 G, carbon molecules are very weakly bound relative to individual atoms). Iron molecules 
are not bound at 75 < 10 13 G, but become energetically more favorable than individual atoms at 
, larger field strengths. 

in ' 

PACS numbers: 31.15.Ew, 95.30.Ky, 97.10.Ld 

> 

\Q ' I. INTRODUCTION 

\o : 

Neutron stars (NSs) are endowed with magnetic fields far beyond the reach of terrestrial laboratories 0, [3, [1]. 
Most of the ~ 1600 known radio pulsars have surface magnetic fields in the range of 10 11 — 10 13 G, as inferred 
\Q \ from their measured spin periods and period derivatives and the assumption that the spindown is due to magnetic 
dipole radiation. A smaller population of older, millisecond pulsars have B ~ 10 8 — 10 9 G. For about a dozen 
accreting neutron stars in binary systems, electron cyclotron features have been detected, implying surface fields of 
B ~ 10 12 — 10 13 G. An important development in astrophysics in the last decade centered on the so-called anomalous 
x-ray pulsars and soft gamma repeaters Q: there has been mounting observational evidence in recent years that 
supports the idea that these are magnetars, neutron stars whose radiations are powered by superstrong magnetic 
> ■ fields of order 10 14 — 10 15 G or higher @, [y, 0|- By contrast, the highest static magnetic field currently produced 
in a terrestrial laboratory is 5 x 10 5 G; transient fields approaching 10 9 G have recently been generated during 
• • ■ high- intensity laser interactions with dense plasmas @. 

It is well-known that the properties of matter can be drastically modified by strong magnetic fields found on 
neutron star surfaces. The natural atomic unit for the magnetic field strength, Bo, is set by equating the electron 
cyclotron energy hu>Be = fr(eB/m e c) = 11.577 B%2 keV, where B12 = B/(10 12 G), to the characteristic atomic energy 
1 e 2 /ao = 2x 13.6 eV (where ao is the Bohr radius): 

B Q = — V~ = 2.3505 x 10 9 G. (1) 

For b = B I Bo ^ 1, the usual perturbative treatment of the magnetic effects on matter (e.g., Zecman splitting of atomic 
energy levels) does not apply. Instead, in the transverse direction (perpendicular to the field) the Coulomb forces act 
as a perturbation to the magnetic forces, and the electrons in an atom settle into the ground Landau level. Because of 
the extreme confinement of the electrons in the transverse direction, the Coulomb force becomes much more effective 
in binding the electrons along the magnetic field direction. The atom attains a cylindrical structure. Moreover, it is 
possible for these elongated atoms to form molecular chains by covalent bonding along the field direction. Interactions 
between the linear chains can then lead to the formation of three-dimensional condensed matter [JJ lid 1 1 11 ] . 

Our main motivation for studying matter in such strong magnetic fields arises from the importance of understanding 
neutron star surface layers, which play a key role in many neutron star processes and observed phenomena. The- 
oretical models of pulsar and magnetar magnetospheres depend on the cohesive properties of the surface matter in 
strong magnetic fields [L?! fl3L IT3. ITEI. liH] . For example, depending on the cohesive energy of the surface matter, an 
acceleration zone ( "polar gap" ) above the polar cap of a pulsar may or may not form. More importantly, the surface 
layer directly mediates the thermal radiations from neutron stars. The advent of x-ray telescopes in recent years 



2 



has made detailed study of neutron star surface emission a reality. Such study can potentially provide invaluable 
information on the physical properties and evolution of NSs: equation of state at supernuclear densities, superfluidity, 
cooling history magnetic field, surface composition, different NS populations, etc. (see, e.g., Ref. More than two 

dozen isolated neutron stars (including radio pulsars, radio-quiet neutron stars and magnetars) have clearly detected 
thermal surface emission 0, [lU, [l9j] . While some neutron stars show featureless spectra, absorption lines or features 
have been detected in half a dozen or so systems Indeed, many of the observed neutron stars have sufficiently low 
surface temperatures and/or high magnetic fields, such that bound atoms or molecules are expected to be present in 
their atmospheres [2CJ, |2l|, |22j, |23[ . It is even possible that the atmosphere is condensed into a solid or liquid form from 
which radiation directly emerges [111, HH, |24| . Thus, in order to properly interpret various observations of neutron 
stars, it is crucial to have a detailed understanding of the properties of atoms, molecules and condensed matter in 
strong magnetic fields (B ~ 10 n -I0 16 G). 

A. Previous works 

H and He atoms at almost all field strengths have been well studied [13, [H|,j26j], including the nontrivial effect 
associated with the center-of-mass motion of a H atom [13] ■ Neuhauser et al. |28l ] presented numerical results for 
several atoms up to Z = 26 (Fe) at B ~ I0 12 G based on calculations using a one-dimensional Hartree-Fock method 
(see also Ref. [29( for Z up to 10). Some results [based on a two-dimensional (2D) mesh Hartree-Fock method] for 
atoms (up to Z = 10) at the field strengths B/Bq = 0.5 — 10 4 are also available [30ll3ll 32]. The Hartree-Fock method 
is approximate because electron correlations are neglected. Due to their mutual repulsion, any pair of electrons tend 
to be more distant from each other than the Hartree-Fock wave function would indicate. In zero-field, this correlation 
effect is especially pronounced for the spin-singlet states of electrons for which the spatial wave function is symmetrical. 
In strong magnetic fields (B 3> Bo), the electron spins (in the ground state) are all aligned antiparallel to the magnetic 
field, and the multielectron spatial wave function is antisymmetric with respect to the interchange of two electrons. 
Thus the error in the Hartree-Fock approach is expected to be less than the 1% accuracy characteristic of zero-field 
Hartree-Fock calculations 12a. [331. Other calculations of heavy atoms in strong magnetic fields include Thomas-Fermi 
type statistical models [34l l35l J36j] and density functional theory [53, [H, [H, Ho] • The Thomas- Fermi type models are 
useful in establishing asymptotic scaling relations, but are not adequate for obtaining accurate binding and excitation 
energies. The density functional theory can potentially give results as accurate as the Hartree-Fock method after 
proper calibration is made [U . 

Quantitative results for the energies of hydrogen molecules Hjv with N = 2, 3, 4, 5 in a wide range of field strengths 
(B ^> Bo) were obtained (based on the Hartree-Fock method) by Lai et al. [Tl|, [43| and molecular excitations were 
studied in Refs. [44],[45j] (more complete references can be found in Ref. [HI). Quantum Monte Carlo calculations of 
H2 in strong magnetic fields have been performed [46j . Some numerical results of He2 for various field strengths are 
also available [ll| . Hartree-Fock results of diatomic molecules (from H2 up to C2) and several larger molecules (up 
to H 5 and He 4 ) at B/B = 1000 are given in Ref. [13]. 

B. Plan of this paper 

In this paper and its companion paper [48j | , we develop a density-functional-theory calculation of the ground-state 
energy of matter for a wide range of magnetic field strengths, from 10 12 G (typical of radio pulsars) to 2 x 10 15 G 
(magnetar fields). We consider H, He, C, and Fe, which represent the most likely composition of the outermost layer 
of neutron stars (e.g., Ref. [3j). The present paper focuses on atoms (and related ions) and small molecules. Because 
of additional complications related to the treatment of band structure, calculations of infinite molecular chains and 
condensed matter are presented in Ref. [48j . 

Our calculations are based on density functional theory [4!|[5(|[IU- As mentioned above, the Hartree-Fock method 
is expected to be highly accurate, particularly in the strong field regime where the electron spins are aligned with each 
other. In this regime the density functional method is not as accurate, due to the lack of an exact correlation function 
for electrons in strong magnetic fields. However, in dealing with systems with many electrons, the Hartree-Fock method 
becomes increasingly impractical as the magnetic field increases, since more and more Landau orbitals (even though 
electrons remain in the ground Landau level) are occupied and keeping track of the direct and exchange interactions 
between electrons in various orbitals becomes computationally rather tedious. Our density-functional calculations 
allow us to obtain the energies of atoms and small molecules and the energy of condensed matter using the same 
method, thus providing reliable cohesive energy of condensed surface of ma gnet ic neutron stars, a main goal of our 
study. Compared to previous density-functional-theory calculations [3?l l38l. l39l. |40| , we use an improved exchange- 
correlation function for highly magnetized electron gases, we calibrate our density functional code with previous 
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results (when available) based on other methods, and (for calculations of condensed matter) adopt a more accurate 
treatment of the band structure. Moreover, our calculations extend to the magnetar-like field regime (B ~ fO 15 G). 

Note that in this paper we neglect the motions of the nuclei, due to electron-nucleus interactions or finite temper- 
atures. The center-of-mass motions of the atoms and molecules induce the motional Stark effect, which can change 
the internal structure of the bound states (see, e.g., Refs. [HI, H3|). Such issues are beyond the scope of this paper. 

After summarizing the approximate scaling relations for atoms and molecules in strong magnetic fields in Sec. II, 
we describe our method in Sec. Ill and present numerical results in Sec. IV. Some technical details are given in the 
Appendix. 



II. BASIC SCALING RELATIONS FOR ATOMS AND MOLECULES IN STRONG MAGNETIC FIELDS 



A. Atoms 



First consider a hydrogenic atom (with one electron and nuclear charge Z). In a strong magnetic field with 
b = B/Bq 3> Z 2 , the electron is confined to the ground Landau level ("adiabatic approximation"), and the Coulomb 
potential can be treated as a perturbation. The energy spectrum is specified by two quantum numbers, (m, v), where 
m = 0, 1, 2, . . . measures the mean transverse separation between the electron and the nucleus (— m is also known as 
the magnetic quantum number), while v specifies the number of nodes in the z wave function. There are two distinct 
types of states in the energy spectrum E mu . The "tightly bound" states have no node in their z wave functions {v = 0). 
The transverse size of the atom in the (m,0) state is L±_ ~ p m = (2m + l) 1 / 2 /^, with po = (Hc/eB) 1 ^ 2 = & -1 / 2 (in 
atomic units). 1 For p m <C 1, the atom is elongated with L z ^> L±. We can estimate the longitudinal size L z by 
minimizing the energy, E ~ L~ 2 — ZL~ X hn(L z /L±) (where the first term is the kinetic energy and the second term 
is the Coulomb energy), giving 

L z ~ (zin-^—) . (2) 

\ Zp m J 

The energy is given by 



E m ~ — Z 



Z 2 V 2m + 1 



(3) 



for b 3> (2m + 1)Z 2 . Another type of state of the atom has nodes in the z wave functions [y > 0). These states are 
"weakly bound", and have energies given by E mv ~ — Z 2 n~ 2 Ry, where n is the integer part of (y + l)/2. The sizes 
of the wave functions are p m perpendicular to the field and L z ~ v 2 jZ along the field (see Ref. and references 
therein for more details). 

A multielectron atom (with the number of electrons N e and the charge of the nucleus Z) can be constructed by 
placing electrons at the lowest available energy levels of a hydrogenic atom. The lowest levels to be filled are the 
tightly bound states with v — 0. When a n /Z ^> ^2N e — lp , i.e., b 3> 2Z 2 N e , all electrons settle into the tightly 
bound levels with m = 0, 1, 2, • • • , N e — 1. The energy of the atom is approximately given by the sum of all the 
eigenvalues of Eq. Accordingly, we obtain an asymptotic expression for N e ^> 1 [52J: 

^-^ e (ln^) 2 . (4) 

For intermediate-strong fields (but still strong enough to ignore Landau excitations), Z 2 N e 2 ' 3 « 6 « 2Z 2 N e , 
many v > states of the inner Landau orbitals (states with relatively small m) are populated by the electrons. In this 
regime a Thomas-Fermi type model for the atom is appropriate, i.e., the electrons can be treated as a one-dimensional 
Fermi gas in a more or less spherical atomic cell (see, e.g., Refs. [Hj, HH). The electrons occupy the ground Landau 
level, with the z momentum up to the Fermi momentum pp ~ n/b, where n is the number density of electrons inside 
the atom (recall that the degeneracy of a Landau level is eB /he ~ b). The kinetic energy of electrons per unit volume 
is £fc ~ bp% ~ n 3 /b 2 , and the total kinetic energy is E k ~ R 3 n 3 /b 2 ~ N 3 /(b 2 R 6 ), where R is the radius of the atom. 



Unless otherwise specified, we use atomic units, in which length is in ao (Bohr radius), mass in m e , energy in e 2 /ao = 2 Ry, and field 
strength in units of Bg. 
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The potential energy is E p ~ —ZN e /R (for N e < Z). Therefore the total energy of the atom can be written as 
E ~ N^/(b 2 R 6 ) - ZN e /R. Minimizing E with respect to R yields 

R~ (N*/Zy/ 5 b- 2 / 5 , E~ -{Z 2 N e f' 5 b 2 / 5 . (5) 

For these relations to be valid, the electrons must stay in the ground Landau level; this requires Z/R <C Tilobc = b, 
which corresponds to b ^> Z 2 N e 2 ^ 3 . 



B. Molecules 

In a strong magnetic field, the mechanism of forming molecules is quite different from the zero-field case 0, [43j. 
Consider hydrogen as an example. The spin of the electron in a H atom is aligned antiparallel to the magnetic field 
(flipping the spin would cost UuiBe): therefore two H atoms in their ground states (m = 0) do not bind together 
according to the exclusion principle. Instead, one H atom has to be excited to the m = 1 state. The two H atoms, one 
in the ground state (m = 0) , another in the m = 1 state then form the ground state of the H2 molecule by covalent 
bonding. Since the activation energy for exciting an electron in the H atom from the Landau orbital m to (m + 1) 
is small, the resulting H2 molecule is stable. Similarly, more atoms can be added to form H3, H4, .... The size of 
the H2 molecule is comparable to that of the H atom. The interatomic separation a and the dissociation energy D 
of the H2 molecule scale approximately as a ~ (Info) -1 and D ~ (ln6) 2 , although D is numerically smaller than the 
ionization energy of the H atom. 

Consider the molecule Zjy, formed out of N neutral atoms Z (each with Z electrons and nuclear charge Z). For 
sufficiently large b (see below), the electrons occupy the Landau orbitals with m = 0, 1, 2, . . . , NZ — 1, and the 
transverse size of the molecule is Lj_ ~ (N Z /b) 1 / 2 . Let a be the atomic spacing and L z ~ Na the size of the molecule 
in the z direction. The energy per "atom" in the molecule, E — En/N, can be written as E ~ Z(Na)~ 2 — (Z 2 /a)l, 
where I ~ ln(a/Lj_). Variation of E with respect to a gives 



a~(ZN 2 l)- L , E~-Z d N 2 l 2 , with I ~ In [j^j ■ (6) 

This above scaling behavior is valid for 1 <C N <C N s . The "critical saturation number" N s is reached when a ~ L±, 
or when W3| 



b 



1/5 



N -[zs) • (7) 

Beyond N s , it becomes energetically more favorable for the electrons to settle into the inner Landau orbitals (with 
smaller m) with nodes in their longitudinal wave functions (i.e., v 7^ 0). For N > N s , the energy per atom asymptotes 

to a value E Z 9 / 5 b 2 / 5 , and size of the atom scales as L± ~ a ~ Z 1 / 5 ^ 2 / 5 , independent of iV — the molecule 

essentially becomes one-dimensional condensed matter. 

The scaling relations derived above are obviously crude — they are expected to be valid only in the asymptotic 
limit, ln(6/Z 3 ) 3> 1. For realistic neutron stars, this limit is not quite reached. Thus these scaling results should 
only serve as a guide to the energies of various molecules. For a given field strength, it is not clear from the above 
analysis whether the Z n molecule is bound relative to individual atoms. To answer this question requires quantitative 
calculations. 



III. DENSITY-FUNCTIONAL CALCULATIONS: METHODS AND EQUATIONS 

Our calculations will be based on the "adiabatic approximation," in which all electrons are assumed to lie in the 
ground Landau level. For atoms or molecules with nucleus charge number Z , this is an excellent approximation for 
b 3> Z 2 . Even under more relaxed condition, b Z 4 ^ 3 (assuming the number of electrons in each atom is N e ~ Z) 
this approximation is expected to yield a reasonable total energy of the system and accurate results for the energy 
difference between different atoms and molecules; a quantitative evaluation of this approximation in this regime is 
beyond the scope of this paper (but see Refs. [30L l3lL l33|). 

In the adiabatic approximation, the one-electron wave function ("orbital") can be separated into a transverse 
(perpendicular to the external magnetic field) component and a longitudinal (along the magnetic field) component: 



= W m (r_i_)/ mv (z) 



(8) 
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Here W m is the ground-state Landau wave function [55| given by 



i ( p /v 

p V2^!VV2po/ XP V 4 Po 



W m (ri_) = 7 == ( -!=— ) exp ( ^ ) exp(-im^) , (9) 



where po — (Hc/eB) 1 / 2 is the cyclotron radius (or magnetic length), and f mv is the longitudinal wave function which 
must be solved numerically. We normalize / m „ over all space: 

dz\f rn „(z)\ 2 = 1, (10) 

J — oc 

so that J dr \^ mi/ (v)\ 2 — 1. The density distribution of electrons in the atom or molecule is 

n(r) = £ |* m „(r)| 2 = ]T \f m u(z)\ 2 \W m \ 2 (p) , (11) 

where the sum is over all the electrons in the atom or molecule, with each electron occupying an (mv) orbital. The 
notation |W m | 2 (p) = | W^, (r ) | 2 is used here because W m is a function of p and 4> but \W m \ 2 is a function only of p. 
In an external magnetic field, the Hamiltonian of a free electron is 

1 / e \ 2 heB , . 

H= — (p + -A) +- a,, (12) 

Zm e V c / 2m e c 

where A = |B x r is the vector potential of the external magnetic field and a z is the z component Pauli spin matrix. 

For electrons in Landau levels, with their spins aligned parallel/antiparallel to the magnetic field, the Hamiltonian 
becomes 

H = ^-+(n L + iy. Be ± 1 -n. Be , (13) 

where nj, = 0, 1,2, • • • is the Landau level index; for electrons in the ground Landau level, with their spins aligned 
antiparallel to the magnetic field (so til — and a z — > —1), 

h = St- ( 14 ) 

2m e 

The total Hamiltonian for the atom or molecule then becomes 

where the sum is over all electrons and V is the total potential energy of the atom or molecule. From this we can 
derive the total energy of the system. 

Note that we use nonrelativistic quantum mechanics in our calculations, even when Tilubc ^ m e c 2 or B > Bq = 
Bq/o 2 = 4.414 x 10 13 G. This is valid for two reasons: (i) The free-electron energy in relativistic theory is 



E = 



c 2 p 2 z +m 2 e c i ( l + 2n L 



B 
~B~a 



"I 1/2 

(16) 



For electrons in the ground Landau level (til — 0), Eq. (fT6|) reduces to E ~ m e c 2 +p 2 /(2m e ) for p z c -C m e c 2 ; the 
electron remains nonrelativistic in the z direction as long as the electron energy is much less than m e c 2 ; (ii) Eq. (J9j) 
indicates that the shape of Landau transverse wave function is independent of particle mass, and thus Eq. ^ is valid 
in the relativistic theory. Our calculations assume that the longitudinal motion of the electron is nonrelativistic. This 
is valid at all field strengths and for all elements considered with the exception of iron at B > 10 15 G. Even at 
B = 2 x 10 15 G (the highest field considered in this paper), however, we find that the most-bound electron in any 
Fe atom or molecule has a longitudinal kinetic energy of only ~ 0.2m e c 2 and only the three most-bound electrons 
have longitudinal kinetic energies > 0.1m e c 2 . Thus relativistic corrections are small in the field strengths considered 
in this paper. Moreover, we expect our results for the relative energies between Fe atoms and molecules to be much 
more accurate than the absolute energies of either the atoms or the molecules. 
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Consider the molecule Zjv, consisting of N atoms, each with an ion of charge Z and Z electrons. In the lowest- 
energy state of the system, the ions are aligned along the magnetic field. The spacing between ions, a, is chosen to 
be constant across the molecule. In the density functional theory, the total energy of the system can be represented 
as a functional of the total electron density n(r): 



E[n] = E K [n] + E eZ [n] + E dil [n] + E cxc [n] + E zz [n] . 



(17) 



Here EK[n] is the kinetic energy of a system of noninteracting electrons, and E eZ , Edit, and E zz are the electron-ion 
Coulomb energy, the direct electron-electron interaction energy, and the ion-ion interaction energy, respectively, 



JV-1 

E zz [n] = Y,(, N -3) 



|r — r 

Z 2 e 2 



J a 



The location of the ions in the above equations is represented by the set {zj}, with 

Zj - = (2j-iV-l)|z. 
The term E oxc represents exchange-correlation energy. In the local approximation, 

E exc [n]= J dr n(r) e cxc (n) , 



(18) 
(19) 
(20) 

(21) 
(22) 



where e cxc (n) = e cx (n) + e corr (n) is the exchange and correlation energy per electron in a uniform electron gas of 
density n. For electrons in the ground Landau level, the (Hartree-Fock) exchange energy can be written as follows 



e C x(™) 



where the dimensionless function Fit) is 



F{t) 



4 / dx 
lo 



tan 



2 p 2 nF(t) , 



In 1 



and 



t = ( ) = 2^pln i 



(23) 



(24) 



(25) 



[tib = (V2k 2 Pq) 1 is the density above which the higher Landau levels start to be filled in a uniform electron gas]. 
For small t, F(t) can be expanded as follows [57]: 



F(t) ~ 3 - 7 - ln4t + | (f - 7 - ln4t) + ^ - 7 - hi At j + 0<? Int), 



(26) 



where 7 = 0.5772 • • • is Euler's constant. We have found that the condition t <C 1 is well satisfied everywhere for 
almost all molecules in our calculations. The notable exceptions are the carbon molecules at B = 10 12 G and the iron 
molecules at B = 10 13 G, which have t < 1 near the center of the molecule. These molecules are expected to have 
higher t values than the other molecules in our calculations, as they have large Z and low B? 



For the uniform gas model, t oc Z 6 / 5 N e B 3 / 5 . 
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The correlation energy of uniform electron gas in strong magnetic fields has not be calculated in general, e xcep t in 
the regime t <C 1 and Fermi wavenumber k F — 2n 2 p 2 n ^> 1 [or n ^> (27r 3 /9ga )^ 1 ]. Skudlarski and Vignale I58J use 
the random-phase approximation to find a numerical fit for the correlation energy in this regime (see also Ref. [591]): 

e corr = -— [0.595(i/6) 1/8 (l - 1.009i 1/8 )] . (27) 
Po 

In the absence of an "exact" correlation energy density we employ this strong-ficld-limit expression. Fortunately, 
because we are concerned mostly with finding energy changes between different states of atoms and molecules, the 
correlation energy term does not have to be exact. The presence or the form of the correlation term has a modest 
effect on the atomic and molecular energies calculated but has very little effect on the energy difference between them 
(see Appendix B for more details on various forms of the correlation energy and comparisons). 

Variation of the total energy with respect to the total electron density, SE[n]/Sn = 0, leads to the Kohn-Sham 
equations: 



2m, 



V 2 + K ff (r) 



# m „(r) = e mv ^ mv {r) , (28) 



where 



N 



3 = 1 

with 



V eS (r) = , _ | + e 2 / dr' — + /w(n), (29) 



d(ne cxc ) 

Mcxc(n) = — q • (30) 



dn 

Averaging the Kohn-Sham equations over the transverse wave function yields a set of one-dimensional equations: 

d2 \W m \ 2 {p) , , U ^ j,\W m \\p)n{v>) 



1m e dz 2 ^-^ 



Ze 



+ / dv ± \W m \\p) («) fmu(z) — £mvfmv\Z) • (^1) 



These equations are solved self-consistently to find the eigenvalue £ m „ and the longitudinal wave function f m u(z) for 
each orbital occupied by the ZN electrons. Once these are known, the total energy of the system can be calculated 
using 



E[n]=J2^mu-^ I I dvdr' n ( r H^) + f drn(r)[e exc (n) - M exc(n)] + ( N ~ ■ ( 32 ) 

mv l r r I J - = l J a 

Details of our method used in computing the various integrals and solving the above equations are given in Appendix 
A. 

Note that for a given system, the occupations of electrons in different (mv) orbitals are not known a priori, and 
must be determined as part of the procedure of finding the minimum energy state of the system. In our calculation, 
we first guess no, n\, ri2, ■ ■ ., the number of electrons in the v = 0, 1, 2, . . . orbitals, respectively (e.g., the electrons in 
the v = orbitals have m — 0, 1, 2, . . . , tiq — 1). Note that uq + n\ + ri2 + ■ ■ ■ = NZ. We find the energy of the system 
for this particular set of electron occupations. We then vary the electron occupations and repeat the calculation until 
the true minimum energy state is found. Obviously, in the case of molecules, we must vary the ion spacing a to 
determine the equilibrium separation and the the ground-state energy of the molecule. Graphical examples of how 
the ground state is chosen are given in the next section. 

IV. RESULTS 

In this section we present our results for the parallel configuration of Hat (up to N = 10), He^v (up to N = 8), 
Cn (up to N = 5), and Fe^r (up to N = 3) at various magnetic field strengths between B = 10 12 G and 2 x 10 15 G. 
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For each molecule (or atom), data is given in tabular form on the molecule's ground-state energy, the equilibrium 
separation of the ions in the molecule, and its orbital structure (electron occupation numbers tiq, ri\, n%, . . .). In 
some cases the first-excited-state energies are given as well, when the ground-state and first-excited-state energies are 
similar in value. We also provide the ground-state energies for selected ionization states of C and Fe atoms; among 
other uses, these quantities are needed for determining the ion emission from a condensed neutron star surface [48| . 
All of the energies presented in this section are calculated to better than 0.1% numerical accuracy (see Appendix A). 

For each of the molecules and ions presented in this section we provide numerical scaling relations for the ground- 
state energy as a function of magnetic field, in the form of a scaling exponent (3 with En oc -Bf 2 - We have provided 
this information to give readers easy access to energy values for fields in between those listed in the tables. The 
ground-state energy is generally not well fit by a constant f3 over the entire magnetic field range covered by this work, 
so we have provided (3 values over several different magnetic field ranges. Note that the theoretical value (3 = 2/5 (see 
Sec. II) is approached only in certain asymptotic limits. 

We discuss here briefly a few trends in the data: All of the molecules listed in the following tables are bound. The 
Fe2 and Fe3 molecules at B\ 2 = 5 are not bound, so we have not listed them here, but we have listed the Fe atom 
at this field strength for comparison with other works. All of the bound molecules listed below have ground-state 
energies per atom that decrease monotonically with increasing N, with the exception of Ha? at B\ 2 = 1, which has a 
slight upward glitch in energy at H4 (see Table QJ. Additionally, these energies approach asymptotic values for large 
N — the molecule essentially becomes one-dimensional condensed matter [481 ] . The equilibrium ion separations also 
approach asymptotic values for large N, but there is no strong trend in the direction of approach: sometimes the 
equilibrium ion separations increase with increasing N, sometimes they decrease, and sometimes they oscillate back 
and forth. 

In general, we find that for a given molecule (e.g., Fe^), the number of electrons in v > states decreases as the 
magnetic field increases. This is because the characteristic transverse size po oc B~ x l 2 decreases, so the electrons 
prefer to stay in the v = states. For a given field strength, as the number of electrons in the system NN e increases 
(e.g., from Fe 2 to Fes), more electrons start to occupy the v > states since the average electron- nucleus separation 
p m ex (2m + l) 1 / 2 ^? -1 / 2 becomes too large for large to. For large enough N the value of no, the number of electrons 
in v = states, levels off, approaching its infinite chain value (see Ref. [481]). Similar trends happen with m, n 2 , etc., 
though much more slowly. 

There are two ways that we have checked the validity of our results by comparison with other works. First, we 
have repeated several of our atomic and molecular calculations using the correlation energy expression empirically 
determined by Jones (37j : 

e 2 

e corr = (0.0096 In pin + 0.122). (33) 

Po 

The results we then obtain for the atomic ground-state energies agree with those of Jones [13, 38]. For example, for 
Fe at B12 = 5 we find an atomic energy of —108.05 eV and Jones gives an energy of —108.18 eV. The molecular 
ground-state energies per atom are of course not the same as those for the infinite chain from Jones's work, but 
they are comparable, particularly for the large molecules. For example, we find for lies a t B X2 = 5 that the energy 
per atom is —1242 eV and Jones finds for He^ that the energy per cell is —1260 eV. (See Appendix B for a brief 
discussion of why in our calculations we chose to use the Skudlarski-Vignale correlation energy expression over that 
of Jones.) 

Second, we have compared our hydrogen, helium, and carbon molecule results to those of Refs. [H, [47} • Because 
these works use the Hartree-Fock method, we cannot compare absolute ground-state energies with theirs, but we can 
compare energy differences. We find fair agreement, though the Hartree-Fock results are consistently smaller. Some 
of these comparisons are presented in the following subsections. 

A. Hydrogen 

Our numerical results for H are given in Table U and Table HT1 Note that at B\ 2 = 1, H4 is less bound than 
H3, and thus E = E^/N is not a necessarily a monotonically decreasing function of N at this field strength. For 
the H4 molecule, two configurations, (no,ri\) = (4,0) and (3, 1), have very similar equilibrium energies (see Fig. [1]), 
although the equilibrium ion separations are different. The real ground state may therefore be a "mixture" of the two 
configurations; such a state would presumably give a lower ground-state energy for H4, and make the energy trend 
monotonia 

Hartree-Fock results for H molecules are given in [43]. For H2, H3, and H4, the energies (per atom) are, respectively: 
-184.3, -188.7, -185.0 eV at B 12 = 1; -383.9, -418.8, -432.9 eV at B 12 = 10; and -729.3, -847.4, -915.0 eV at 
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TABLE I: Ground-state energies, ion separations, and electron configurations of hydrogen molecules, over a range of magnetic 
field strengths. In some cases the first-excited-state energies are also listed. Energies are given in units of eV, separations in 
units of ao (the Bohr radius). For molecules (Hjv) the energy per atom is given, E — En/N. All of the H and H2 molecules 
listed here have electrons only in the v — states. For the H3 and larger molecules here, however, the molecular structure is 
more complicated, and is designated by the notation (no, ni, . . .), where no is the number of electrons in the v = orbitals, ni 
is the number of electrons in the u = 1 orbitals, etc. 



B12 


H 

E 


H 2 
E a 


H 3 

E a (n ,ni) 


H 4 

E a (n ,ni) 


H 5 

E a (n ,ni) 


1 


-161.4 


-201.1 0.25 


-209.4 0.22 (3,0) 
-191.1 0.34 (2,1) 


-208.4 0.21 (4,0) 
-207.9 0.26 (3,1) 


-213.8 0.23 (4,1) 
-203.1 0.200 (5,0) 


10 


-309.5 


-425.8 0.125 


-469.0 0.106 (3,0) 


-488.1 0.096 (4,0) 


-493.5 0.090 (5,0) 
-478.9 0.112 (4,1) 


100 


-540.3 


-829.5 0.071 


-961.2 0.057 (3,0) 


-1044.5 0.049 (4,0) 


-1095.5 0.044 (5,0) 


1000 


-869.6 


-1540.5 0.044 


-1818.0 0.033 (3,0) 


-2049 0.028 (4,0) 


-2222 0.024 (5,0) 



B\2 


H 6 

E a (no,ni) 


H 8 

E a (no,ni,n2) 


Hio 

E a (no,ni,n2) 


1 


-214.1 0.23 (4,2) 
-213.4 0.21 (5,1) 


-215.8 0.23 (5,2,1) 
-215.3 0.25 (4,3,1) 


-216.2 0.22 (6,3,1) 
-216.0 0.23 (5,3,2) 


10 


-496.5 0.101 (5,1) 
-490.8 0.86 (6,0) 


-507.1 0.095 (8,2,0) 
-504.1 0.089 (7,1,0) 


-509.3 0.091 (7,3,0) 
-506.8 0.087 (8,2,0) 


100 


-1125.0 0.041 (6,0) 


-1143.0 0.038 (8,0,0) 
-1139.5 0.043 (7,1,0) 


-1169.5 0.038 (9,1,0) 
-1164.0 0.042 (8,2,0) 


1000 


-2351 0.22 (6,0) 


-2518 0.0190 (8,0,0) 


-2600 0.0170 (10,0,0) 
-2542 0.0200 (9,1,0) 



TABLE II: Fit of the ground-state energies of hydrogen molecules to the scaling relation E oc Bf 2 . The scaling exponent [3 is 
fit for each molecule Hjv over three magnetic field ranges: B\2 = 1 — 10, 10 — 100, and 100 — 1000. 



B12 




H H2 H3 H4 H5 He Hg Hio 


1-10 


0.283 0.326 0.350 0.370 0.363 0.365 0.371 0.372 


10-100 


0.242 0.290 0.312 0.330 0.346 0.355 0.353 0.361 


100-1000 


0.207 0.269 0.277 0.293 0.307 0.320 0.343 0.347 



B12 = 100. Thus, our density-functional-theory calculation tends to overestimate the energy \E\ by about 10%. Note 
that the Hartree-Fock results also reveal a non-monotonic behavior of E at N — 4 for B\ 2 = 1, in agreement with 
our density-functional result. Demeur et al. [47| calculated the energies of H2-H5 at B12 = 2.35; their results exhibit 
similar trends. 



B. Helium 

Our numerical results for He are given in Table IIIII and Table IIVI 

The energies (per atom) of He and He2 based on Hartree-Fock calculations [ll[ are, respectively: —575.5, —601.2 eV 
at B 12 = 1; -1178, -1364 eV at B 12 = 10; -2193, -2799 eV at B 12 = 100; and -3742, -5021 eV at B 12 = 1000. 
At B\ 2 = 2.35, Demeur et al. [47] find that the energies (per atom) of He, He2, He3, and He4 are, respectively: 
— 753.4, —812.6, —796.1, —805.1 eV. Using our scaling relations, we find for that same field that the energies of He, 
He2, He3, and Hes (we do not have an He4 result) are: —791, —871, —889, —901 eV. Thus, our density-functional 
theory calculation tends to overestimate the energy \E\ by about 10%. 
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FIG. 1: Molecular energy per atom versus ion separation for various hydrogen molecules at B12 = 1. The energy of the H atom 
is shown as a horizontal line at —161.4 eV. The two lowest-energy configurations of H4 have nearly the same minimum energy, 
so the curves for both configurations are shown here. 



C. Carbon 



Our numerical results for C are given in Table [V] Table fVTI and Table Ivnl 

The only previous result of C molecules is that by Demeur et al. [47[, who calculated C2 only at Bn = 2.35. At 
this field strength, our calculation shows that C2 is bound relative to C atom [E = —5994, —6017 eV for C, C2), 
whereas Demeur et al. find no binding (E — —5770, —5749 eV for C, C2). Thus our result differs qualitatively from 
[47I ]. We also disagree on the ground-state occupation at this field strength: we find (no,ni) = (9,3) while Demeur 
et al. find (rig, «i) = (7, 5). We suggest that if Demeur et al. used the occupation (no, ni) — (9, 3) they would obtain 
a lower-energy for C2, though whether C2 would then be bound remains uncertain. Since the numerical accuracy of 
our computation is 0.1% of the total energy (thus, about 6 eV for B12 = 2.35), our results for B12 < a few should be 
treated with caution. 

Figure [5] gives some examples of the longitudinal electron wave functions. One wave function of each node type 
in the molecule (y = to 4) is represented. Note that on the atomic scale each wave function is nodeless in nature; 
that is, there are no nodes at the ions, only in between ions. The exception to this is at the central ion, where due to 
symmetry considerations the antisymmetric wave functions must have nodes. [The nodes for (m, v) — (0, 2) are near, 
but not at, the ions j ' — 2 and j = 4. This is incidental.] This is not surprising when one considers that all of the 
electrons in atomic carbon at this field strength are nodeless. The entire molecular wave function can be thought of as 
a string of atomic wave functions, one around each ion, each modified by some phase factor to give the overall nodal 
nature of the wave function. Indeed, for atoms at field strengths that are low enough to allow v > states, we find 
that their corresponding molecules have electron wave functions with nodes at the ions. Atomic Fe at B\2 = 10, for 
example, has an electron wave function with one node at the ion, and Fe2 at B\2 = 10 has an electron wave function 
with a node at each ion. 
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TABLE III: Ground-state energies, ion separations, and electron configurations of helium molecules, over a range of magnetic 
field strengths. In some cases the first-excited-state energies are also listed. Energies are given in units of eV, separations in 
units of ao (the Bohr radius). For molecules (Hejv) the energy per atom is given, E = En/N. All of the He and He2 molecules 
listed here have electrons only in the v — states. For the He3 and larger molecules here, however, the molecular structure is 
more complicated, and is designated by the notation (no, m, . . •), where no is the number of electrons in the v = orbitals, ni 
is the number of electrons in the v = 1 orbitals, etc. 



B\2 


ne 
E 


xie2 
E a 


ne3 

E a (no,ni) 


ne5 

E a (no,ni,n 2 ) 


nes 

E a (n , m, n 2 ,n 3 ) 


1 


-603.5 


-641.2 0.25 


-647.3 0.28 (5,1) 
-633.0 0.32 (4,2) 


-653.1 0.29 (6,3,1) 
-649.4 0.28 (7,2,1) 


-656.7 0.28 (7,5,3,1) 
-656.5 0.27 (8,5,2,1) 


10 


-1252.0 


-1462.0 0.115 


-1520.0 0.105 (6,0) 
-1462.0 0.125 (5,1) 


-1553.5 0.110 (8,2,0) 
-1547.5 0.105 (9,1,0) 


-1574.5 0.110 (10,5,1,0) 
-1574.0 0.105 (11,4,1,0) 


100 


-2385 


-3039 0.060 


-3370 0.050 (6,0) 
-3140 0.054 (5,1) 


-3573 0.044 (10,0,0) 
-3543 0.049 (9,1,0) 


-3694 0.045 (13,3,0,0) 
-3690 0.043 (14,2,0,0) 


1000 


-4222 


-5787 0.036 


-6803 0.028 (6,0) 


-7887 0.022 (10,0,0) 


-8406 0.0200 (15,1,0,0) 
-8357 0.0180 (16,0,0,0) 



TABLE IV: Fit of the ground-state energies of helium molecules to the scaling relation E <x B^ 2 ■ The scaling exponent f5 is fit 
for each molecule Hejv over three magnetic field ranges: B\2 = 1 — 10, 10 — 100, and 100 — 1000. 



B\2 


P 

He He 2 He3 Hes Hes 


1-10 


0.317 0.358 0.371 0.376 0.380 


10-100 


0.280 0.318 0.346 0.362 0.370 


100-1000 


0.248 0.280 0.305 0.344 0.357 



D. Iron 

Our numerical results for Fe are given in Table IVIII| Table IIX1 and Table [X] The energy curves for Byi = 500 are 
shown in Fig. El and some results for B\i = 100 are shown in Fig. |H 

There is no previous quantitative calculation of Fe molecules in strong magnetic fields that we are aware of. The 
most relevant work is that of Abrahams and Shapiro [60] , who use a Thomas- Fermi type model to calculate Fe and 
Fe2 energies for magnetic fields up to B\2 = 30. Unfortunately, a comparison of our results with those of this work is 
not very useful, as Thomas-Fermi models are known to give inaccurate energies and in particular large overestimates 



TABLE V: Ground-state energies, ion separations, and electron configurations of carbon molecules, over a range of magnetic 
field strengths. In some cases the first-excited-state energies are also listed. Energies are given in units of eV, separations in 
units of ao (the Bohr radius). For molecules (Cjv) the energy per atom is given, E = En/N. All of the C atoms listed here 
have electrons only in the v = orbitals. For the C2 and larger molecules here, however, the molecular structure is more 
complicated, and is designated by the notation (no, ni, . . .), where no is the number of electrons in the v — orbitals, m is the 
number of electrons in the v = 1 orbitals, etc. 





C 


c 2 






c 3 






c 4 






c 5 


B\2 


E 


E a 


(n ,ni) 


E 


a 


(n ,ni,n 2 ) 


E 


a 


(no, rn, 712,113) 


E 


a 


(no, ni, 112,713,714) 


1 


-4341 


-4351 0.53 


(8,4) 


-4356 


0.52 


(9,6,3) 


-4356 


0.52 


(10,7,4,3) 


-4358 


0.48 


(11,8,6,3,2) 






-4349 0.46 


(9,3) 


-4354 0.50 


(10,5,3) 


-4354 0.56 


(9,7,5,3) 


-4357 0.47 


(12,8,5,3,2) 


10 


-10075 


-10215 0.150 


(11,1) 


-10255 


0.175 


(13,4,1) 


-10255 


0.180 


(15,6,2,1) 


-10275 


0.150 


(18,8,3,1) 






-10200 0.180 


(10,2) 


-10240 


0.185 


(14,3,1) 


-10250 


0.185 


(14,7,2,1) 


-10270 


0.155 


(17,9,3,1) 


100 


-21360 


-23550 0.054 


(12,0) 


-24060 


0.055 


(17,1,0) 


-24350 


0.054 


(21,3,0,0) 


-24470 


0.057 


(23,6,1,0,0) 










-23960 


0.058 


(16,2,0) 


-24300 


0.056 


(20,4,0,0) 


-24460 


0.056 


(24,5,1,0,0) 


1000 


-41330 


-50760 0.027 


(12,0) 


-54870 


0.024 


(18,0,0) 


-56500 


0.024 


(23,1,0,0) 


-57640 


0.022 


(28,2,0,0,0) 
















-56190 


0.022 


(24,0,0,0) 


-57520 


0.023 


(27,3,0,0,0) 
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TABLE VI: Ground-state energies of ionized carbon atoms over a range of magnetic field strengths. Energies are given in units 
of eV. For these field strengths, the electron configuration of C atoms is such that all of their electrons lie in the v = orbitals; 
therefore the ionized atoms have all electrons in the v = orbitals as well. The ionization state is designated by the notation, 
"C n+ ," where n is the number of electrons that have been removed from the atom. The entry "C 5+ ," for example, is a carbon 
nucleus plus one electron. 



Bl2 


C 


C+ 


C 2+ 


C 3+ 


C 4+ 


C 5+ 


1 


-4341 


-4167 


-3868 


-3411 


-2739 


-1738.0 


10 


-10075 


-9644 


-8917 


-7814 


-6213 


-3877 


100 


-21360 


-20370 


-18730 


-16300 


-12815 


-7851 


1000 


-41330 


-39210 


-35830 


-30920 


-24040 


-14425 



TABLE VII: Fit of the ground-state energies of neutral and ionized carbon atoms and carbon molecules to the scaling relation 
E oc i?f 2 . The scaling exponent f3 is fit over three magnetic field ranges: B12 = 1 — 10, 10 — 100, and 100 — 1000. 



B12 


C 5+ C 4+ c+ 


P 
C 


C2 C3 C4 C5 


1-10 


0.348 0.356 0.364 


0.366 


0.371 0.372 0.372 0.372 


10-100 


0.306 0.314 0.325 


0.326 


0.363 0.370 0.376 0.377 


100-1000 


0.264 0.273 0.284 


0.287 


0.334 0.358 0.366 0.372 




-0.3 ^ 1 1 1 1 1 L - 

10 20 30 

z(Po) 



FIG. 2: Longitudinal wave functions for selected electron orbitals of C5 at B12 = 1, at the equilibrium ion separation. Different 
orbitals are labeled by (m,u). Only the 2 > region is shown. Wave Functions with even v are symmetric about 2 = 0, and 
those with odd v are antisymmetric about 2 = 0. The filled circles denote the ion locations. 



13 



TABLE VIII: Ground-state energies, ion separations, and electron configurations of iron molecules, over a range of magnetic 
field strengths. In some cases the first-excited-state energies are also listed. Energies are given in units of keV, separations in 
units of ao (the Bohr radius). For molecules (Fejv) the energy per atom is given, E — En/N. The electron configuration is 
designated by the notation (no, ni, . . .), where no is the number of electrons in the v = orbitals, m is the number of electrons 
in the v = 1 orbitals, etc. Note that no information is listed for the Fe2 and Fe3 molecules at B12 = 5, as we have found that 
these molecules are not bound at this field strength. Also note that there are two lowest-energy states for Fe2 at B12 = 100; 
within the error of our calculation, the two states have the same minimum eneriges. 





Fe 




Fe 2 




Fe 3 




B\2 


E 


(no, ni) 


E 


a 


(no,ni) 


E a 


(no, ni, n2) 


5 


-107.20 


(24,2) 






10 


-142.15 


(25,1) 


-142.18 


0.30 


(32,19,1) 




100 


-354.0 


(26,0) 


-354.9 


0.107 


(39,13) 


-355.2 0.107 


(47,21,10) 








-354.9 


0.103 


(40,12) 


-355.1 0.108 


(46,22,10) 


500 


-637.8 


(26,0) 


-645.7 


0.048 


(45,7) 


-648.1 0.048 


(58,16,4) 








-645.4 


0.050 


(44,8) 


-648.0 0.050 


(57,16,5) 


1000 


-810.6 


(26,0) 


-828.8 


0.035 


(47,5) 


-834.1 0.035 


(62,13,3) 








-828.4 


0.034 


(48,4) 


-834.0 0.036 


(61,14,3) 


2000 


-1021.5 


(26,0) 


-1061.0 


0.025 


(49,3) 


-1073.0 0.025 


(67,10,1) 








-1056.0 


0.023 


(50,2) 


-1072.5 0.025 


(66,11,1) 



TABLE IX: Ground-state energies of ionized iron atoms over a range of magnetic field strengths. Energies are given in units of 
keV. For B12 > 100, the electron configuration of Fe atoms is such that all of their electrons lie in the v — orbitals; therefore 
for these field strengths the ionized atoms have all electrons in the v = orbitals as well. The ionization state is designated 
by the notation, "Fe n+ ," where n is the number of electrons that have been removed from the atom. The entry "Fe 25+ ," for 
example, is an iron nucleus plus one electron. 



B\2 


Fe 


Fe+ 


Fe 2+ 


Fe 3 + 


Fe 4 + 


Fe 5 + 


Fe 10 + 


Fe 15+ 


Fe 20+ 


Fe 25+ 


100 


-354.0 


-352.8 


-351.2 


-349.0 


-346.4 


-343.2 


-318.3 


-273.8 


-199.65 


-59.01 


500 


-637.8 


-635.3 


-632.0 


-627.8 


-622.7 


-616.6 


-569.4 


-486.5 


-350.2 


-99.48 


1000 


-810.6 


-807.2 


-802.8 


-797.2 


-790.7 


-782.5 


-715.8 


-602.0 


-439.6 


-122.70 


2000 


-1021.5 


-1016.0 


-1008.5 


-999.8 


-989.1 


-976.7 


-905.4 


-768.6 


-546.8 


-150.10 



of binding and cohesive energies. As an example, from Ref. [601 ] the energy difference between Fe and Fe2 at Byi = 30 
is 1.7 keV, which is twice as large as our result at B\2 — 100. 

In Table |VlTT] we have not provided results for the Fe2 and Fe3 molecules at B\2 = 5, as these molecules are not 
bound relative to the Fe atom. We have not provided results for the Fe3 molecule at B12 — 10 because the energy 
difference (per atom) between Fe3 and the Fe atom at this field strength is smaller than the error in our calculation, 
0.1% of \E\ or 140 eV. The energy difference (per atom) between the Fe2 molecule and the Fe atom at B12 = 10 is 
also smaller than the error in our calculation (indeed, the difference should be less than that between Fe3 and Fe at 
this field strength), but we have redone the calculation using more grid and integration points such that the energy 
values reported here for these two molecules are accurate numerically to 0.01% (see Appendix A). At this accuracy, 
our results indicate that Fe2 is bound over Fe at B12 = 10 with a energy difference per atom of 30 eV. 

Figure [4] illustrates how the ground-state electron configuration is found for each molecule. The configuration with 
the lowest equilibrium energy is chosen as the ground-state configuration. In the case depicted in Fig. [U Fe2 at 
B\2 = 100, there are actually two such configurations. Within the error of our calculation, we cannot say which 
one represents the ground state. Note that the systematic error seen in the minimization curves of the various Fe2 
configurations is much smaller than our target 0.1% error for the total energy (the sinusoidal error in the figure has 
an amplitude of » 30 eV, or around 0.01% of the total energy). 

V. CONCLUSIONS 

We have presented density-functional-theory calculations of the ground-state energies of various atoms and molecular 
chains (Hjv up to N = 10, Hejv up to N = 8, Cjv up to N = 5, and Fe^r up to N = 3) in strong magnetic fields ranging 
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TABLE X: Fit of the ground-state energies of neutral and ionized iron atoms and iron molecules to the scaling relation E oc B 
The scaling exponent is fit over three magnetic field ranges: B12 = 100 — 500, 500 — 1000, and 1000 — 2000. 



B12 


Fe 25+ p e 20+ Fe 10+ p e + 



Fe 


Fe 2 Fe 3 


100-500 


0.324 0.349 0.361 0.365 


0.366 


0.372 0.374 


500-1000 


0.303 0.328 0.330 0.345 


0.346 


0.359 0.364 


1000-2000 


0.291 0.315 0.339 0.332 


0.334 


0.358 0.363 




from B = 10 12 G to 2 x 10 15 G. These atoms and molecules may be present in the surface layers of magnetized neutron 
stars, such as radio pulsars and magnetars. While previous results (based on Hartree-Fock or density-functional-theory 
calculations) are available for some small molecules at selected field strengths (e.g., Refs. [ll|> Hj, \41§) many other 
systems (e.g., larger C molecules and Fe molecules) are also computed in this paper. We have made an effort to present 
our numerical results systematically, including fitting formulae for the independence of the energies. Comparison with 
previous results (when available) show that our density-functional calculations tend to overestimate the binding energy 
l-Ejvl by about 10%. Since it is advantageous to use the density functional theory to study systems containing large 
number of electrons (e.g., condensed matter; see Ref. |4a|). it would be useful to find ways to improve upon this 
accuracy. 

At B12 > 1, hydrogen, helium, and carbon molecules are all more energetically favorable than their atomic coun- 
terparts (although for carbon, the relative binding between the atom and molecule is rather small), but iron is not. 
Iron molecules start to become bound at B12 ^ 10, and are not decidedly more favorable than isolated atoms until 
about B 12 = 100. 

For the bound molecules considered here, the ground-state energy per atom approaches an asymptotic value as N 
gets large. The molecule then essentially becomes a one-dimensional infinite chain. We will study such condensed 
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-354 



-354.2 
| -354.4 
\ -354.6 

LU 

-354.8 



-355 

0.1 0.105 0.11 0.115 

a(a ) 

FIG. 4: Molecular energy per atom versus ion separation for various configurations of electrons in the Fe2 molecule at B12 = 100. 
The configurations are labeled using the notation (no, ni), where no is the number of electrons with v = and ni is the number 
with v = 1. The energy of the Fe atom is shown as a horizontal line at —354.0 keV. The states "(40, 12)" and "(39, 13)" have 
the lowest equilibrium energies of all possible configurations and within the numerical accuracy of our calculation have the 
same equilibrium energies. The wavy structure of the curves gives an indication of the numerical accuracy of our code. Note 
that states with electrons in the v — 2 orbitals [for example, (39, 12, 1)] have energies higher than the atomic energy and are 
therefore unbound. 

matter in our companion paper [48| . 
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APPENDIX A: NUMERICAL METHOD 

1. Evaluating the integrals in the Kohn-Sham equations 

The two most computation-intensive terms in the Kohn-Sham equations [Eq. (|31[) ] are the ion-electron interaction 
term and the direct electron-electron interaction term: 




(Al) 
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and 



^ee,m (^) 



dr± dr' 



, \W m \ 2 (p)n(v') 



(A2) 



Equation (|A1[) . together with the exchange-correlation term, J dr± \W m \ 2 (p) j!/ exc (n), can be integrated by a standard 
quadrature algorithm, such as Romberg integration [61] . Equation (|A2|) . however, is more complicated and its 
evaluation is the rate-limiting step in the entire energy calculation. The integral is over four variables (p, p' , z' , and 
cf> or (/) — 0'), so it requires some simplification to become tractable. To simplify the integral we use the identity (see, 
e.g., Ref. [621 ]) 



OO /> oo 

V / dqe m (^J n ( q p)J n ( qp >)e-^- z '\ 



where J n (z) is the nth order Bessel function of the first kind. Then 

j <v') 



V ee (r) = J dr' 



r - r' 



/oo p oo r poo 

dz' I dqJ (qp) / p' dp' n(p' , z')J (qp r ) 
-oo Jo L^o 



exp(-gr|z - z'\) 



and 



Ke,m(2) = / dr ± \W m \ 2 (p)V ee (r) 



/oo /»oo r /-c 

dz' dq 
-oo JO UO 



pdp\W m \ 2 {p) J (qp) 



p dp n(p',z')J Q (qp') 



exp(— q\z — z'\ 



Using Eq. (fTTj) for the electron density distribution, Eq. (IA5I) becomes 

/oo />oo 
d^'|/mv(2')| 2 / G m (q)G m < (q) exp(-g|z - z'\) , 
-oo JO 



where 



and 



G m {q) = 2tt / pdp\W m \ 2 (p)J (qp) 
JO 

= exp(- <Z 2 /2)i m (<Z 2 /2), 



M*) = -73-sr(* m O 



to! dx r 

is the Laguerre polynomial of order m. These polynomials can be calculated using the recurrence relation 

mL m (x) = (2to - 1 - x)L m ^i(x) - (to - l)L m _ 2 (x) , 



(A3) 



(A4) 



(A5) 



(A6) 



(A7) 



(A8) 



(A9) 



with Lo(x) = 1 and Li(x) = 1 — x. 

Using the method outlined above the original four-dimensional integral in Eq. (|A2[) reduces to a two-dimensional 
integral. Once a value for z is specified, the integral can be evaluated using a quadrature algorithm (such as the 
Romberg integration method). 



2. Solving the differential equations and total energy 



The Kohn-Sham equations [Eq. (|3"T|)] are solved on a grid in z. Because of symmetry we only need to consider 
z > 0, with z = at the center of the molecule. The number and spacing of the z grid points determine how 
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accurately the equations can be solved. In this paper we have attempted to calculate ground-state energies to better 
than 0.1% numerical accuracy. This requires approximately (depending on Z and B) 133 grid points for a single atom 
calculation, plus 66 more for each additional atom in the molecule, or a total of w 66 * (N + 1) points for an iV-atom 
molecule. The grid spacing is chosen to be constant from the center out to the outermost ion, then exponentially 
increasing as the potential decays to zero. The maximum z value for the grid is chosen such that the amplitude of all 
of the electron wave functions / m „ at that point is less than 5 x 1CP 3 . 

For integration with respect to p, p', or q (e.g., when calculating the direct electron-electron interaction term), our 
0.1%-accuracy goal for the energy values requires an accuracy of approximately 10~ 5 in the integral. A variable-step- 
size integration routine is used for each such integral, where the number of points in the integration grid is increased 
until the error in the integration is within the desired accuracy. 

Solving the Kohn-Sham equations requires two boundary conditions for each (mi') orbital. The first is that f mv (z) 
vanishes exponentially for large z. Because the j m v(z) wave functions must be symmetric or anti-symmetric about 
the center of the molecule, there is a second boundary condition: wave functions with an even number of nodes 
have an extremum at the center and wave functions with an odd number of nodes have a node at the center; i.e., 
fm V (0) = for even v and / m! /(0) = for odd v. In practice, we integrate Eq. (|31|) from the large- z edge of the z 
grid and "shoot" toward z = 0, adjusting e mv until the boundary condition at the center is satisfied. One final step 
must be taken to ensure that we have obtained the desired energy and wave function shape, which is to count the 
number of nodes in the wave function. For each (mis) orbital there is only one wave function shape that satisfies the 
required boundary conditions and has the correct number of nodes v (e.g., the shape of each orbital in Fig. [21 however 
complicated- looking, is uniquely determined). 

To determine the electronic structure of an atom or molecule self-consistently, a trial set of wave functions is first 
used to calculate the potential as a function of z, and that potential is used to calculate a new set of wave functions. 
These new wave functions are then used to find a new potential, and the process is repeated until consistency is 
reached. In practice, we find that f mv (z) = works well as the trial wave function, and rapid convergence can be 
achieved: four or five iterations for atoms and no more than 20 iterations for the largest and most complex molecules. 
To prevent overcorrection from one iteration to the next, the actual potential used for each iteration is a combination 
of the newly-generated potential and the old potential from the previous iteration (the weighting used is roughly 30% 
old, 70% new). 

APPENDIX B: CORRELATION ENERGY 

As was mentioned in Sec. Ill, the form of the correlation energy has very little effect on the relative energy between 
atom and molecule (or between different states of the same molecule). This holds true even if the calculations are 
done in the extreme case where the correlation energy term is set to zero. As an example, consider the energy of the 
C2 molecule at B = 10 15 G. Using the correlation energy of Skudlarski and Vignale [Eq. (|27)) ]. we find the C atom 
has energy E a — —41 330 eV and the C2 molecule has energy per atom E m = —50 760 eV, so that the relative energy 
is AE = 9430 eV. Using the correlation energy of Jones [Eq. (|3"3"])]. we find E a = -44 420 eV and E m = -53 840 eV, 
so that AE = 9420 eV. Without any correlation term at all, E a = -38 600 eV and E m = -47 960 eV, so that 
AE = 9360 eV. As another example of the relative unimportance of the correlation term, two other works using 
density- functional calculations, Jones [37j . Relovsky and Ruder [40| . find very similar cohesive energy (i.e., infinite 
chain) results even though they use two very different correlation energy terms. For example, at B = 5 x 10 12 G they 
both find a cohesive energy of 220 eV for Heoo. 

We make one final comment about the accuracy of our chosen correlation energy term, the Skudlarski- Vignale 
expression Eq. (|2"T)) . Jones [13] found an empirical expression for the correlation energy at high B [Eq. (|3"3")) ]. and then 
checked its accuracy using the fact that the self- interaction of an occupied, self-consistent orbital should be zero, i.e., 

Edir[n m v] + ficxcKJ = , (Bl) 

where n mu = |4' m!/ (r)| 2 is the number density of electrons in the (mv) orbital. Performing such a test on the 
Skudlarski- Vignale expression, we find that the error in Eq. (|B1[) is of order 5-20% for B\2 — 1 and up to 20-30% 
for B\2 = 1000 for the elements and molecules considered here. Testing Jones's expression, we find it does as well 
and sometimes better at Bn = 1, but at large fields it does considerably worse, up to 60-100% error for B12 = 1000. 
For example, for He2 at B12 — 1000 the Skudlarski- Vignale correlation function satisfies Eq. (|B1[) to within 23% but 
Jones's expression satisfies Eq. (|B1[) only to within 63%. Thus, the Skudlarski- Vignale correlation function adopted 
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in this paper is much more accurate than Jones's expression for a wide range of field strengths. 
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